clear; clc; close all;
T = 300;
k0 = 235.988436 - 0.663191454*T^1 + 0.00403072646*T^2 - 1.39161256e-5*T^3 + 2.39071664e-8*T^4 - 1.98206006e-11*T^5 + 6.36310053e-15*T^6;
fprintf('计算结果 k0 = %.8f (T = %d K)\n', k0, T);       
kf = 0.026;    
ks = k0;         
c = linspace(0, 1, 1001);
keff = zeros(size(c));  

for i = 1:length(c)
    ci = c(i); 
    if 0 <= ci && ci <= 0.15
        numerator = ks * (kf + 2*ks + 2*ci*(kf - ks));
        denominator = kf + 2*ks - ci*(kf - ks);  
        keff(i) = numerator / denominator;
    elseif 0.15 < ci && ci <= 0.65
        A = kf*(3*ci - 1) + ks*(2 - 3*ci);
        keff(i) = 0.25 * (A + sqrt(A^2 + 8*ks*kf));
    elseif ci > 0.65

        term1 = 1/(kf - ks);   
        term2 = (1 - ci)/(3*ks);
        keff(i) = ks + ci / (term1 + term2);
    end
end

figure('Color','white','Position',[100,100,800,600]);
plot(c, keff, 'LineWidth',2, 'Color','#2E86AB');      

hold on;
plot([0.15,0.15], [min(keff),max(keff)], '--r', 'LineWidth',1.5); 
plot([0.65,0.65], [min(keff),max(keff)], '--r', 'LineWidth',1.5);  
scatter([0.15,0.65], [keff(find(c==0.15)), keff(find(c==0.65))], ...
        60, 'red', 'filled');  
title(['有效导热系数随体积分数的变化曲线 (T = ', num2str(T), ' K)'], 'FontSize',16, 'FontName','SimHei'); 

xlabel('体积分数 c', 'FontSize',14, 'FontName','SimHei');  
ylabel('有效导热系数 k_{eff} (W/(m·K))', 'FontSize',14, 'FontName','SimHei'); 

legend(['k_{eff}-c曲线 (T=', num2str(T), ' K)'], '分段点', 'Location','best', 'FontSize',12, 'FontName','SimHei');  

text(0.8, max(keff)*0.9, ['温度: ', num2str(T), ' K'], 'FontSize',14, 'FontName','SimHei', 'Color','#E63946', 'FontWeight','bold');

grid on; 
grid minor;
set(gca, 'FontSize',12, 'FontName','SimHei');  
ylim([min(keff)*0.95, max(keff)*1.05]); 

fprintf('c=0.15时，keff=%.4f (T = %d K)\n', keff(find(c==0.15)), T);
fprintf('c=0.65时，keff=%.4f (T = %d K)\n', keff(find(c==0.65)), T);